Introduction

https://github.com/edwardlavender/patter

Introduction

Team

Institutions

Workshop

• The goal today is to introduce patter
• We will focus on how to use the package
• There will be:
  • Interactive presentation & questions
  • (Time allowing) Optional practical activity (using example datasets or your own)

Context

• patter uses a dedicated high-performance backend written in Julia

https://github.com/edwardlavender/Patter.jl

Evolution

• patter evolved from flapper (Lavender et al., 2023)

Evolution

• flapper is a more general-purpose package for passive acoustic telemetry that supports:
  • data acquisition, simulation, assembly and processing
  • spatial operations and distance metrics (e.g., fast C++ LCP algorithms)
  • data exploration (detection statistics, movement metrics)
  • ‘flapper’ algorithms

Evolution

• patter focuses specifically on model inference for SSMs using particle algorithms

https://github.com/edwardlavender/patter

Operating systems

• patter v.1.0.0 supports Windows and MacOS
• patter@dev supports Linux and will be rolled out soon

Installation

• To install patter, follow the instructions on https://github.com/edwardlavender/patter
  • Install R and Julia
  • On Windows, install RTools
  • On MacOS, install Xcode (and configure .Makevars as required)
  • Install patter (and let it handle Julia dependencies)

Connect to Julia

• The syntax for julia_connect() is:

julia_connect(..., JULIA_PROJ, .threads, ...)

• In patter@dev, this is:

julia_connect(JULIA_HOME, JULIA_PROJ, JULIA_NUM_THREADS, ...)

• If you’ve set environmental variables (recommended), just use:

julia_connect()

Overview

## Overview

Overview

• We will work though:
  • Boilerplate setup
  • An example workflow with real data
  • An example simulation-based workflow

Boilerplate setup

• Connect to Julia:

julia_connect()

• Set the seed in R and Julia via set_seed():

set_seed(2024 - 10 - 18)

Study system

Observations include detections at receivers:

head(dat_acoustics)

   individual_id           timestamp receiver_id
           <int>              <POSc>       <int>
1:            25 2016-03-17 01:50:00          26
2:            25 2016-03-17 01:52:00          26
3:            25 2016-03-17 01:54:00          26
4:            25 2016-03-17 01:58:00          26
5:            25 2016-03-17 02:00:00          26
6:            25 2016-03-17 02:04:00          26

Study system

We also have archival (depth) observations, collected at a resolution of two-minutes:

head(dat_archival)

   individual_id           timestamp  depth
           <num>              <POSc>  <num>
1:            25 2016-03-16 00:00:00 159.15
2:            25 2016-03-16 00:02:00 161.01
3:            25 2016-03-16 00:04:00 162.86
4:            25 2016-03-16 00:06:00 164.48
5:            25 2016-03-16 00:08:00 158.46
6:            25 2016-03-16 00:10:00 155.45

Study system

• Our goal in the first part of this presentation is to reconstruct movements and patterns of space use for an example individual given these data.

# Define detections
det <-
  dat_acoustics |>
  filter(individual_id == 25L) |>
  mutate(individual_id = NULL) |>
  as.data.table()

# Define archival (depth) observations
arc <-
  dat_archival |>
  filter(individual_id == 25L) |>
  mutate(individual_id = NULL) |>
  as.data.table()

Study system

plot(arc$timestamp, arc$depth * -1, 
     col = "royalblue", type = "l",
     ylim = c(-max(arc$depth), 0), 
     xlab = "Time stamp", ylab = "Depth (m)")
points(det$timestamp, rep(0, nrow(det)))

Study system

• We need to define a map of our study system in R and Julia

Study system

Study system

• Cell entries should be floats
• Code inhospitable habitats (e.g., land) as NA_real_
• A planar (UTM) coordinate reference system is required

Study system

• Define the study timeline:
  • This sets the period, and resolution, at which we represent animal movement
  • We may have observations regularly or irregularly along this timeline

Study system

• Use assemble_timeline() to set a timeline based on the input datasets:

timeline <- assemble_timeline(.datasets = list(det, arc), 
                              .step = "2 mins", 
                              .trim = TRUE)
str(timeline)

 POSIXct[1:24225], format: "2016-03-17 01:50:00" "2016-03-17 01:52:00" "2016-03-17 01:54:00" ...

Study system

• Set the timeline using seq.POSIXt():

timeline <- seq(as.POSIXct("2016-03-17 01:50:00", tz = "UTC"),
                as.POSIXct("2016-03-18 01:48:00", tz = "UTC"), 
                by = "2 mins")
str(timeline)

 POSIXct[1:720], format: "2016-03-17 01:50:00" "2016-03-17 01:52:00" "2016-03-17 01:54:00" ...

State

• State is an Abstract Type in Patter.jl
• This is an object that holds the parameters describing the ‘state’ of an individual
• Typically, state means ‘location’ (in two or three dimensions), but additional parameters may be included

# The Julia code for the StateXY struct
struct StateXY <: State
    map_value::Float64
    x::Float64
    y::Float64
end 

State

• In Patter.jl, the following sub-types are built-in:

  • StateXY(map_value, x, y): Used for two dimensional (x, y) states ;
  • StateXYZD(map_value, x, y, z, angle): Used for four-dimensional (x, y, z, direction) states;

State

• In patter, you just need to specify the state sub-type as a character string:
  • "StateXY" maps onto StateXY in Julia
  • "StateXYZD" maps onto StateXYDZ in Julia

state <- "StateXY"

• Each state-type maps onto a corresponding movement model

Movement model

• ModelMove is an Abstract Type in Patter.jl
• This is an object that holds the dimensions of the movement model, e.g.,
  • The map
  • A distribution of step lengths
  • A distribution of turning angles

# The Julia code for the ModelMoveXY struct
struct ModelMoveXY{T, U, V} <: ModelMove
    map::T
    dbn_length::U
    dbn_angle::V
end 

Movement model

• The movement model can be instantiated via an R move_*() wrapper
  • StateXY maps onto ModelMoveXY, which is wrapped by move_xy()
  • StateXYZD maps onto ModelMoveXYZD, which is wrapped by move_xyzd()
• The fields of the movement model in the move_*() wrapper must be specified using Julia code

Movement model

• Define a two-dimensional random walk (x, y):

# Define maximum moveable speed between two time steps
mobility   <- 750.0

# Define 2D RW (XY)
model_move <- move_xy(dbn_length = glue("truncated(Gamma(1.0, 250.0), upper = {mobility})"),
                      dbn_angle = "Uniform(-pi, pi)")

Movement model

• move_*() wrappers return a string of Julia code:

model_move

ModelMoveXY(env, truncated(Gamma(1.0, 250.0), upper = 750), Uniform(-pi, pi));

Movement model

• Visualise the components of our movement model:

expr <- expression({
  pp <- par(mfrow = c(1, 2))
  curve(dtrunc(x, spec = "gamma", a = 0, b = mobility, shape = 1.0, scale = 250.0), 
        from = 0, to = mobility, n = 1e3L,
        xlab = "Step length (m)", ylab = "Density")
  curve(dunif(x, -pi, pi),
        from = -pi - 0.1, to = pi + 0.1, n = 1e3L,
        xlab = "Heading (rad)", ylab = "Density")
  par(pp)
})

Movement model

• Visualise the components of our movement model:

Movement model

• Visualise realisations of the movement model:

Movement model

• More sophisticated movement models are possible, including:
  • Random walks in two or three dimensions
  • Correlated random walks
  • Behavioural-switching models
• See the help for ?State to get started

Movement model

• This is an example of a more sophisticated movement model for flapper skate:

Example behavioural-switching model.

Observation model(s)

• We have defined the movement model
• This encodes the process of animal movement
• We now need to define the observation models
• These encode the probability of the observations, given the latent (unobserved) location of the animal

Observation model(s)

• As a reminder, we have collected detections and depth time series:

head(det, 2)

             timestamp receiver_id
                <POSc>       <int>
1: 2016-03-17 01:50:00          26
2: 2016-03-17 01:52:00          26

head(arc, 2)

             timestamp  depth
                <POSc>  <num>
1: 2016-03-16 00:00:00 159.15
2: 2016-03-16 00:02:00 161.01

Observation model(s)

• Each dataset should comprise a time series of observations for one individual
• Datasets should be pre-processed as required: e.g.,
  • Check receiver locations are valid on the map
  • Identify false detections
• patter accepts any kind of observational dataset

Observation model(s)

• Each dataset should be provided as a data.table with the following columns:
  • sensor_id, timestamp, obs;
  • Additional columns that define observation-model parameters;

Observation model(s)

• In passive acoustic telemetry datasets, we record detections but the observations comprise detections and non-detections
• Use the assemble_acoustics() helper to assemble the acoustic (detection, non-detection) time series:

acc <- assemble_acoustics(.timeline  = timeline,
                          .acoustics = det,
                          .moorings  = dat_moorings, 
                          .services  = NULL)

Observation model(s)

• If you have passive acoustic telemetry observations, you should also use acoustic containers
• The syntax is:

# ModelObsAcousticContainer
containers <- assemble_acoustics_containers(.acoustics, .mobility, ...)

Observation model(s)

• We have to assemble other datasets with the required columns (sensor_id, timestamp, obs, ...) manually

head(arc, 2)

             timestamp  depth
                <POSc>  <num>
1: 2016-03-16 00:00:00 159.15
2: 2016-03-16 00:02:00 161.01

arc <- 
  arc |> 
  # Define 'dummy' sensor_id
  mutate(sensor_id = 1L) |> 
  # Observations must be in the 'obs' column
  rename(obs = "depth") |> 
  # Select columns
  select("sensor_id", "timestamp", "obs") |>
  # Add observation model parameters 
  mutate(depth_sigma = 50, depth_deep_eps = 20) |>
  as.data.table()

Observation model(s)

• It is a good idea to pass other datasets through the assemble_archival() function for processing/checks though

arc <- assemble_archival(.timeline = timeline, .archival = arc)
head(arc)

   sensor_id           timestamp   obs depth_sigma depth_deep_eps
       <int>              <POSc> <num>       <num>          <num>
1:         1 2016-03-17 01:50:00 73.78          50             20
2:         1 2016-03-17 01:52:00 73.32          50             20
3:         1 2016-03-17 01:54:00 73.32          50             20
4:         1 2016-03-17 01:56:00 73.32          50             20
5:         1 2016-03-17 01:58:00 73.55          50             20
6:         1 2016-03-17 02:00:00 68.70          50             20

Observation model(s)

• For each dataset, we need to define the corresponding ModelObs sub-type
• A ModelObs sub-type is a structure that holds the parameters of an observation model
• ModelObs sub-types are instantiated in Julia using the data.table(s) assembled above
• We use ModelObs instances to simulate observations or to evaluate the log-probability of observations, given latent location

Observation model(s)

• The following ModelObs structures are built-in to Patter.jl:
  • ModelObsAcousticLogisTrunc
  • ModelObsDepthNormalTrunc
  • ModelObsDepthUniform

Observation model(s)

• In a truncated logistic detection probability model:
  • receiver_alpha and receiver_beta affect the rate of decline in detection probability with distance from a receiver
  • receiver_gamma is the maximum detection range

Observation model(s)

Example truncated logistic detection probability models

Observation model(s)

• You can explore different parameters values using this code:

expr <- expression({
  # Define parameters 
  receiver_alpha <- 4
  receiver_beta  <- -0.01
  receiver_gamma <- 1000
  # Visualise model 
  curve(ifelse(x <= receiver_gamma, 
              1 / (1 + exp(-(receiver_alpha + receiver_beta * x))), 
              0), 
       from = 0, to = receiver_gamma, 
       xlab = "Distance (m)", ylab = "Detection probability")
})

Observation model(s)

eval(expr)

Observation model(s)

• ModelObsDepthNormalTrunc is a ModelObs structure for a depth observation and a truncated normal model
• This is a suitable depth-observation model for benthic/demersal species
• The individual must be located in an envelope around the bathymetric depth, defined by a normal distribution centred at this location, truncated at depth_deep_eps m below the seabed and 0 m, with standard deviation depth_sigma

Observation model(s)

Example truncated logistic detection probability models

Observation model(s)

• ModelObsDepthUniform is ModelObs structure for a depth observation and a uniform depth model
• This model assumes that an individual must be located in an envelope around the bathymetric depth, defined by two error terms (depth_shallow_eps and depth_deep_eps)
• This model can be tailored for benthic or pelagic species

Observation model(s)

• For the particle filter, multiple sets of observations should be provided as a named list
• List names define the model sub-types

yobs <- list(ModelObsAcousticLogisTrunc = acc, 
             ModelObsDepthNormalTrunc = arc)

Observation model(s)

• As for the movement model, the parameters of the observation model(s) are assumed known

Particle filter

• We have now defined:
  1. Study area map
  2. Study timeline
  3. State
  4. Movement model
  5. Observations/observation model parameters/observation model sub-types
• We are now in a position to run the particle filter!

Two-filter smoother

• Run the two-filter smoother using outputs from pf_filter() in the Julia environment:

smo <- pf_smoother_two_filter(.n_particle = 100L, .n_sim = 100L)

Two-filter smoother

• pf_smoother_two_filter() returns the familiar pf_particles-class object:

class(smo)

[1] "list"         "pf_particles"

summary(smo)

            Length Class      Mode   
xinit       0      -none-     NULL   
states      6      data.table list   
diagnostics 4      data.table list   
convergence 1      -none-     logical

Two-filter smoother

• We can analyse the outputs in the same way as those from the filter:

head(smo$states)

   path_id timestep           timestamp map_value        x       y
     <int>    <int>              <POSc>     <num>    <num>   <num>
1:       1        1 2016-03-17 01:50:00  68.53316 709262.9 6253335
2:       1        2 2016-03-17 01:52:00  75.20494 709235.3 6253358
3:       1        3 2016-03-17 01:54:00  56.53801 709310.5 6253250
4:       1        4 2016-03-17 01:56:00  80.03239 709132.5 6253351
5:       1        5 2016-03-17 01:58:00  84.47292 709281.4 6253504
6:       1        6 2016-03-17 02:00:00  80.03239 709163.0 6253326

Two-filter smoother

• diagnostics includes ESS:

head(smo$diagnostics)

   timestep           timestamp      ess maxlp
      <int>              <POSc>    <num> <num>
1:        1 2016-03-17 01:50:00      NaN   NaN
2:        2 2016-03-17 01:52:00 34.94540   NaN
3:        3 2016-03-17 01:54:00 80.75311   NaN
4:        4 2016-03-17 01:56:00 51.64853   NaN
5:        5 2016-03-17 01:58:00 42.20092   NaN
6:        6 2016-03-17 02:00:00 79.73982   NaN

Two-filter smoother

• convergence is necessarily TRUE:

smo$convergence

[1] TRUE

Mapping

• We can compute probability-of-use (POU) using particle samples from the filter (for comparison) or the smoother (recommended):

map_pou(.map = map, .coord = smo$states)$ud

Mapping

Mapping

• In practice, POU is sensitive to grid resolution and kernel smoothing may be desirable

map_dens(
  .map,
  .owin = as.owin.SpatRaster(.map),
  .coord = NULL,
  .discretise = FALSE,
  .shortcut = list(),
  ...,
  .fterra = FALSE,
  .plot = TRUE,
  .use_tryCatch = TRUE,
  .verbose = getOption("patter.verbose")
)

Mapping

Residency

• We can easily compute residency in a particular area (at any arbitrary time-scale of our choosing) as the proportion of particles inside/outside an area:

# Define an example polygon (e.g., MPA)
spatpoly <- 
  cbind(708433.9, 6251765) |> 
  terra::vect(crs = terra::crs(map)) |> 
  terra::buffer(width = 2000)

# Visualise study area and 'MPA'
terra::plot(map)
points(smo$states$x, smo$states$y, pch = ".", col = "red")
terra::lines(spatpoly, col = "black", lwd = 2)

Residency

Residency

# Lookup whether or not particles are in the polygon
spatcoord <- terra::vect(cbind(smo$states$x, smo$states$y))
smo$states[, in_poly := 
               terra::relate(spatcoord, spatpoly, relation = "intersects")]
               
# Calculate the expected time spent in the 'MPA' (%)
table(smo$states$in_poly) / nrow(smo$states) * 100

   FALSE     TRUE 
71.76111 28.23889 

Boilerplate setup (1/7)

# Load & attach package 
library(patter)

# Load datasets
head(det, 2)
head(arc, 2)

# Connect to Julia 
julia_connect()
set_seed()

# Set map 
map <- dat_gebco()
set_map()

# Study timeline 
timeline <- assemble_timeline(.datasets = list(det, arc))

Set state (2/7)

state <- "StateXY"

Set movement model (3/7)

mobility   <- 750.0
model_move <- move_xy(dbn_length = glue("truncated(Gamma(1.0, 250.0), upper = {mobility})"),
                      dbn_angle = "Uniform(-pi, pi)")

Set observation model(s) (4/7)

# Acoustic observations 
acc <- assemble_acoustics(.timeline  = timeline,
                          .acoustics = det,
                          .moorings  = dat_moorings, 
                          .services  = NULL)

# Archival observations 
arc <- assemble_archival(.timeline = timeline, 
                         .archival = 
                           arc |> 
                             rename(obs = "depth") |> 
                             select("sensor_id", "timestamp", "obs") |>
                             mutate(depth_sigma = 50, depth_deep_eps = 20) |> 
                           as.data.table())

# All observations
yobs <- list(ModelObsAcousticLogisTrunc = acc, 
             ModelObsDepthNormalTrunc = arc)

Run particle filter (5/7)

# List filter arguments
args <- list(.map         = map,
             .timeline    = timeline,
             .state       = state,
             .xinit_pars  = list(mobility = mobility),
             .yobs        = yobs,
             .model_obs   = names(yobs),
             .model_move  = model_move,
             .n_particle  = 1e5L)

# Forward run
fwd <- do.call(pf_filter, args, quote = TRUE)

# Backward run
args$.direction <- "backward"
bwd <- do.call(pf_filter, args, quote = TRUE)

Run particle smoother (6/7)

smo <- pf_smoother_two_filter()

Analyse outputs (7/7)

map_dens(.map = map, .coord = smo, sigma = bw.diggle)

Simulation-based workflow

• patter also supports simulation-based workflows:
  1. sim_array() simulates arrays
  2. sim_path_walk() simulates movement paths
  3. sim_observations() simulates observations
• You can model simulated observations

Simulate arrays

• Use sim_array() to simulate randomly or regularly arranged receivers
• Locations are sampled from .map
• Default observation model parameters (receiver_alpha, receiver_beta, receiver_gamma) can be assigned

moorings <- sim_array(.map = map, 
                      .timeline = timeline,
                      .arrangement = "random",
                      .n_receiver = 100L, 
                      .receiver_alpha = 5, 
                      .receiver_beta = -0.01, 
                      .receiver_gamma = 1200, 
                      .n_array = 1L, 
                      .plot = TRUE)

Simulate arrays

Simulate arrays

head(moorings)

   array_id receiver_id      receiver_start        receiver_end receiver_x
      <int>       <int>              <POSc>              <POSc>      <num>
1:        1           1 2016-03-17 01:50:00 2016-03-18 01:48:00   702042.1
2:        1           2 2016-03-17 01:50:00 2016-03-18 01:48:00   709342.1
3:        1           3 2016-03-17 01:50:00 2016-03-18 01:48:00   707842.1
4:        1           4 2016-03-17 01:50:00 2016-03-18 01:48:00   711042.1
5:        1           5 2016-03-17 01:50:00 2016-03-18 01:48:00   705842.1
6:        1           6 2016-03-17 01:50:00 2016-03-18 01:48:00   708842.1
   receiver_y receiver_alpha receiver_beta receiver_gamma
        <num>          <num>         <num>          <num>
1:    6249607              5         -0.01           1200
2:    6268407              5         -0.01           1200
3:    6269107              5         -0.01           1200
4:    6267107              5         -0.01           1200
5:    6252607              5         -0.01           1200
6:    6253707              5         -0.01           1200

Simulate paths

• We have already seen how to simulate paths

# Set state and movement model 
state      <- "StateXY"
mobility   <- 500.0
model_move <- 
  move_xy(dbn_length = "truncated(Normal(50, 250.0), lower = 0, upper = {mobility})", 
                      dbn_angle = "VonMises(0, 0.25)")

# Simulate path using movement model 
pp <- par(mfrow = c(1, 2))
path_coord <- sim_path_walk(.map = map, 
                            .timeline = timeline, 
                            .state = state, 
                            .n_path = 1L, 
                            .model_move = model_move)

# Fit UD using cross validation
path_ud <- map_dens(.map = map, .coord = path_coord, .discretise = TRUE)$ud
par(pp)

Simulate paths

Simulate paths

head(path_coord)

   path_id timestep           timestamp map_value        x       y
     <int>    <int>              <POSc>     <num>    <num>   <num>
1:       1        1 2016-03-17 01:50:00  54.62202 708042.1 6251707
2:       1        2 2016-03-17 01:52:00  54.17428 708113.0 6251676
3:       1        3 2016-03-17 01:54:00  50.58004 708358.7 6251598
4:       1        4 2016-03-17 01:56:00  48.12802 708078.8 6251506
5:       1        5 2016-03-17 01:58:00  53.05067 708200.9 6251690
6:       1        6 2016-03-17 02:00:00  50.70027 708233.6 6251503

Simulate observations

• To simulate observations, we need to choose our observation models and their parameters
• For each observation, we need to provide parameters as a data.table that defines, for each sensor_id (e.g., receiver), the observation-model parameters associated with that sensor

Simulate observations

• Let’s set things up to simulate acoustic observations:

# ModelObsAcousticLogisTrunc 
# > We have multiple acoustic sensors
# > This data.table includes the relevant parameters
pars_ModelObsAcousticLogisTrunc <-
    moorings |>
    select(sensor_id = "receiver_id",
           "receiver_x", "receiver_y",
           "receiver_alpha", "receiver_beta", "receiver_gamma") |>
    as.data.table()

Simulate observations

• And depth observations:

# ModelObsDepthNormalTrunc
# > We have one depth sensory 
# > This data.table includes the relevant parameters 
pars_ModelObsDepthNormalTrunc <- data.frame(sensor_id = 1L,
                                            depth_sigma = 20,
                                            depth_deep_eps = 20)

Simulate observations

• Simulate observations along .timeline via sim_observations():

obs <- sim_observations(
  .timeline = timeline, 
  .model_obs = c("ModelObsAcousticLogisTrunc",
                 "ModelObsDepthNormalTrunc"),
  .model_obs_pars = list(pars_ModelObsAcousticLogisTrunc, 
                         pars_ModelObsDepthNormalTrunc))

Simulate observations

• sim_observations() returns a list, with one element for every .model_obs
• Each element is a list, with one sub-element for each simulated path
• Each sub-element is a data.table that contains the observations

Simulate observations

str(obs)

List of 2
 $ ModelObsAcousticLogisTrunc:List of 1
  ..$ :Classes 'data.table' and 'data.frame':   72000 obs. of  8 variables:
  .. ..$ timestamp     : POSIXct[1:72000], format: "2016-03-17 01:50:00" "2016-03-17 01:50:00" ...
  .. ..$ obs           : int [1:72000] 0 0 0 0 0 0 0 0 0 0 ...
  .. ..$ sensor_id     : int [1:72000] 1 2 3 4 5 6 7 8 9 10 ...
  .. ..$ receiver_x    : num [1:72000] 702042 709342 707842 711042 705842 ...
  .. ..$ receiver_y    : num [1:72000] 6249607 6268407 6269107 6267107 6252607 ...
  .. ..$ receiver_alpha: num [1:72000] 5 5 5 5 5 5 5 5 5 5 ...
  .. ..$ receiver_beta : num [1:72000] -0.01 -0.01 -0.01 -0.01 -0.01 -0.01 -0.01 -0.01 -0.01 -0.01 ...
  .. ..$ receiver_gamma: num [1:72000] 1200 1200 1200 1200 1200 1200 1200 1200 1200 1200 ...
  .. ..- attr(*, ".internal.selfref")=<externalptr> 
 $ ModelObsDepthNormalTrunc  :List of 1
  ..$ :Classes 'data.table' and 'data.frame':   720 obs. of  5 variables:
  .. ..$ timestamp     : POSIXct[1:720], format: "2016-03-17 01:50:00" "2016-03-17 01:52:00" ...
  .. ..$ obs           : num [1:720] 46.4 14.8 28.5 42.6 61.1 ...
  .. ..$ sensor_id     : int [1:720] 1 1 1 1 1 1 1 1 1 1 ...
  .. ..$ depth_sigma   : num [1:720] 20 20 20 20 20 20 20 20 20 20 ...
  .. ..$ depth_deep_eps: num [1:720] 20 20 20 20 20 20 20 20 20 20 ...
  .. ..- attr(*, ".internal.selfref")=<externalptr> 

Simulate observations

acc <- obs$ModelObsAcousticLogisTrunc[[1]]
head(acc)

             timestamp   obs sensor_id receiver_x receiver_y receiver_alpha
                <POSc> <int>     <int>      <num>      <num>          <num>
1: 2016-03-17 01:50:00     0         1   702042.1    6249607              5
2: 2016-03-17 01:50:00     0         2   709342.1    6268407              5
3: 2016-03-17 01:50:00     0         3   707842.1    6269107              5
4: 2016-03-17 01:50:00     0         4   711042.1    6267107              5
5: 2016-03-17 01:50:00     0         5   705842.1    6252607              5
6: 2016-03-17 01:50:00     0         6   708842.1    6253707              5
   receiver_beta receiver_gamma
           <num>          <num>
1:         -0.01           1200
2:         -0.01           1200
3:         -0.01           1200
4:         -0.01           1200
5:         -0.01           1200
6:         -0.01           1200

Simulate observations

arc <- obs$ModelObsDepthNormalTrunc[[1]]
head(arc)

             timestamp      obs sensor_id depth_sigma depth_deep_eps
                <POSc>    <num>     <int>       <num>          <num>
1: 2016-03-17 01:50:00 46.42963         1          20             20
2: 2016-03-17 01:52:00 14.84555         1          20             20
3: 2016-03-17 01:54:00 28.52349         1          20             20
4: 2016-03-17 01:56:00 42.62775         1          20             20
5: 2016-03-17 01:58:00 61.07769         1          20             20
6: 2016-03-17 02:00:00 34.44895         1          20             20

Simulate observations

# Plot depths
  plot(arc$timestamp, arc$obs * -1, 
       col = "royalblue", type = "l",
       ylim = c(-max(arc$obs), 0), 
       xlab = "Time stamp", ylab = "Depth (m)")
# Add detections 
det <- acc[obs == 1L, ]
points(det$timestamp, rep(0, nrow(det)), col = "red")

Simulate observations

Mapping

pp <- par(mfrow = c(1, 3))

# Plot the UD for the simulated path 
terra::plot(path_ud)
map_hr_home(path_ud, .add = TRUE)

# Plot the estimated UD 
patter_ud <- map_dens(.map = map, 
                   .coord = smo$states, 
                   .discretise = TRUE, 
                   .fterra = TRUE,
                   sigma = bw.diggle)$ud
map_hr_home(patter_ud, .add = TRUE)

# Compare to a UD based on the COA algorithm
coa_ud <- map_dens(.map = map, 
                   .coord = coa(.map = map, 
                                .acoustics = det, 
                                .moorings = moorings, 
                                .delta_t = "2 hours", 
                                .plot_weights = FALSE), 
                   .discretise = FALSE, 
                   .fterra = TRUE, 
                   sigma = bw.diggle)$ud

par(pp)

Mapping

Residency

Parallelisation

Custom States, ModelObs and ModelMove

Custom state -> custom movement model -> wrapper R move_*() function StateXYZ -> ModelMoveXYZ -> move_xyz() states_init.StateXYZ

Patter.simulate_step(state::StateXYZ, model::ModelMoveXYZ, …) logpdf_step(state_from::StateXYZ, state_to::StateXYZ, model_move::ModelMoveXYZ, …)

simulate_obs(state::StateXYZD, model::ModelObsDepthNormal, …) logpdf_obs(state::State, model::ModelObsDepthNormal, …)